#include "../DWMRICore/VectorField.h"
#include "../DWMRICore/Fiber.h"

#include <omp.h>

#include <string>

using namespace std;

char *pVectorFieldPathname = NULL;
char *pStreamlinePathname = NULL;

float stepSize = 0.1f;
float maxLength = 5.0f;
uint innerSteps = 1;

CVectorField *pVectorField = NULL;
uint w = 0;
uint h = 0;
uint d = 0;

void PrintUsage()
{
	printf("Usage: ");
	printf("Streamline -vectorfield <vector_field> -streamline <streamline> [-stepsize <stepsize> -maxlength <maxlength>\n");
	printf("\n");
}

int ParseArguments(const int argc, char *argv[])
{
	for (int i = 1; i < argc; ++i) {
		if (strcmp(argv[i], "-streamline") == 0) {
			pStreamlinePathname = argv[++i];
		} else if (strcmp(argv[i], "-vectorfield") == 0) {
			pVectorFieldPathname = argv[++i];
		} else if (strcmp(argv[i], "-stepsize") == 0) {
			stepSize = (float)atof(argv[++i]);
		} else if (strcmp(argv[i], "-maxlength") == 0) {
			maxLength = (float)atof(argv[++i]);
		} else {
			printf("unknown %s\n", argv[i]);
			return -1;
		}
	}

	if (pVectorFieldPathname == NULL || pStreamlinePathname == NULL || maxLength < 0.0001f || stepSize < 0.00001f)
		return -2;

	return 0;
}

void StreamlineTracking(Seed seed, CFiber *streamline)
{
	seed.dir = pVectorField->GetVector(seed.pos.m_x, seed.pos.m_y, seed.pos.m_z);

	uint maxSteps = uint(maxLength / (stepSize * innerSteps));

	CPoint3F currPos, nextPos;
	CVector3F dir;

	// forward
	currPos = seed.pos;
	streamline->m_pF[0] = currPos;
	streamline->m_fCount = 1;
	dir = seed.dir;
	dir.Normalize();
	for (uint step = 1; step < maxSteps; ++step) {
		bool flag = true;
		for (uint i = 0; i < innerSteps; ++i) {
			nextPos = currPos + dir * stepSize;
			if (nextPos.m_x < 0 || nextPos.m_x > (w-1) ||
				nextPos.m_x < 0 || nextPos.m_x > (w-1) ||
				nextPos.m_x < 0 || nextPos.m_x > (w-1) ) {
					flag = false;
					break;
			}
			currPos = nextPos;
			dir = pVectorField->GetVector(currPos.m_x, currPos.m_y, currPos.m_z);
			dir.Normalize();
		}
		if (flag == false)
			break;

		streamline->m_pF[step] = currPos;
		streamline->m_fCount++;
	}

	// backward
	currPos = seed.pos;
	streamline->m_pB[0] = currPos;
	streamline->m_bCount = 1;
	dir = seed.dir;
	dir.Normalize();
	for (uint step = 1; step < maxSteps; ++step) {
		bool flag = true;
		for (uint i = 0; i < innerSteps; ++i) {
			nextPos = currPos - dir * stepSize;
			if (nextPos.m_x < 0 || nextPos.m_x > (w-1) ||
				nextPos.m_x < 0 || nextPos.m_x > (w-1) ||
				nextPos.m_x < 0 || nextPos.m_x > (w-1) ) {
					flag = false;
					break;
			}
			currPos = nextPos;
			dir = pVectorField->GetVector(currPos.m_x, currPos.m_y, currPos.m_z);
			dir.Normalize();
		}
		if (flag == false)
			break;

		streamline->m_pB[step] = currPos;
		streamline->m_bCount++;
	}
}


int main(int argc, char *argv[])
{
	if (ParseArguments(argc, argv) != 0) {
		PrintUsage();
		return 0;
	}

	pVectorField = new CVectorField();
	pVectorField->OpenVectorFieldFile(pVectorFieldPathname);

	w = pVectorField->GetWidth();
	h = pVectorField->GetHeight();
	d = pVectorField->GetDepth();

	FILE *fp = fopen(pStreamlinePathname, "w");
	fprintf(fp, "%d %d %d %d %f %f\n", w, h, d, 1, stepSize, maxLength);

	std::string pathname = pStreamlinePathname;
	size_t pos = pathname.find_last_of(SLASH);
	std::string path = pathname.substr(0, pos+1);

	for (uint z = 100; z < 150; ++z) {
		printf("Streamline tracking on slice z = %d ... ", z);

		CFiber *streamlines = new CFiber[w*h];
		for (uint y = 0; y < h; ++y) { 
			for (uint x = 0; x < w; ++x) {
				Seed seed;
				seed.pos = CPoint3F((float)x, (float)y, (float)z);

				uint index = y * w + x;
				StreamlineTracking(seed, &(streamlines[index]));
			}
			float ff = float(y+1) / float(h) * 100.0f;
			printf("\rStreamline tracking on slice z = %d ... %.2f%%", z, ff);
		}

		char confPathname[PATHNAME_SIZE];
		char vertPathname[PATHNAME_SIZE];
		sprintf(confPathname, "%sslice_z=%08.04f.conf", path.c_str(), (float)z/(float)1);
		sprintf(vertPathname, "%sslice_z=%08.04f.vert", path.c_str(), (float)z/(float)1);
		fprintf(fp, "%s\t%s\n", confPathname, vertPathname);

		SaveFibers(confPathname, vertPathname, streamlines, w, h);

		SafeDeleteArray(streamlines);

		printf("\rStreamline tracking on slice z = %d ... 100.00%% Done.\n", z);
	}

	SafeDelete(pVectorField);

	return 0;
}